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Abstract 

Three symbolic algorithms for testing the integrability of polynomial systems of 
partial differential and differential-difference equations are presented. The first al- 
gorithm is the well-known Painleve test, which is applicable to polynomial systems 
of ordinary and partial differential equations. The second and third algorithms allow 
one to explicitly compute polynomial conserved densities and higher-order symme- 
tries of nonlinear evolution and lattice equations. 

The first algorithm is implemented in the symbolic syntax of both Macsyma and 
Mathematica. The second and third algorithms are available in Mathematica. The 
codes can be used for computer-aided integrability testing of nonlinear differential 
and lattice equations as they occur in various branches of the sciences and engineer- 
ing. Applied to systems with parameters, the codes can determine the conditions 
on the parameters so that the systems pass the Painleve test, or admit a sequence 
of conserved densities or higher-order symmetries. 
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1 Introduction 



During the last three decades, the study of integrability, invariants, symme- 
tries, and exact solutions of nonlinear ordinary and partial differential equa- 
tions (ODEs, PDEs) and differential- difference equations (DDEs) has been 
the topic of major research projects in the dynamical systems and the soliton 
communities. 

Various techniques have been developed to determine whether or not PDEs 
belong to the privileged class of completely integrable equations [15]. One 
of the most successful and widely applied techniques is the Painleve test, 
named after the French mathematician Paul Painleve (1863-1933) [35], who 
classified second-order differential equations that are globally integrable in 
terms of elementary functions by quadratures or by linearization. In essence, 
the Painleve test verifies whether or not solutions of differential equations in 
the complex plane are single- valued in the neighborhood of all their movable 
singularities. 

To a large extent, the Painleve test is algorithmic, yet very cumbersome when 
done by hand. In particular, the verification of the compatibility conditions is 
assessed by many practitioners of the Painleve test as a painstaking compu- 
tation. In addition to the tedious verification of self-consistency (or compati- 
bility) conditions, computer programs are helpful at exploring all possibilities 
of balancing singular terms. Indeed, the omission of one or more choices of 
"dominant behavior" can lead to wrong conclusions [12,44]. 

We therefore developed the programs painsing. max and painsys.max [25] (both 
in Macsyma syntax [36]), and painsing. m and painsys.m in Mathematica lan- 
guage [53], that perform the Painleve test for polynomial systems of ODEs 
and PDEs. In this paper we demonstrate the above codes by analyzing a 
few prototypical nonlinear equations and systems, such as the Boussinesq and 
nonlinear Schrodinger equations, a class of fifth-order KdV equations, and the 
Hirota-Satsuma and Lorenz systems. 

Our computer code does not deal with the theoretical shortcomings of the 
Painleve test as identified by Kruskal and others [30-32]. Thus far, we have 
implemented the traditional Painleve test [26,27], and not yet incorporated 
the latest advances in Painleve type methods, such as the poly-Painleve test 
[31] or other generalizations [30,32]. Neither did we code the weak Painleve 
test [43,44] or other variants [9-11,17]. Furthermore, we do not have code for 
the singularity confinement method [23], i.e. an adaptation of the Painleve 
test that allows one to test the integrability of difference equations. 

Among the various alternatives to establish the integrability [15] of nonlinear 
PDEs and DDEs, the search for conserved densities and higher-order symme- 



2 



tries is particularly appealing. Indeed, in this paper we will give algorithms 
that apply to both the continuous and semi-discrete cases. We implemented 
these algorithms [18-20] in Mathematica, but they are fairly simple to code in 
other computer algebra languages (see [18,20]). 

Our algorithms are based on the concept of dilation (scaling) invariance. That 
inherently limits their scope to polynomial conserved densities and higher- 
order symmetries of polynomial systems. Although the existence of a sequence 
of conserved densities predicts integrability, the nonexistence of polynomial 
conserved quantities does not preclude integrability. Indeed, integrable PDEs 
or DDEs could be disguised with a coordinate transformation in DDEs that 
no longer admit conserved densities of polynomial type [46]. The same care 
should be taken in drawing conclusions about non-integrability based on the 
lack of higher-order symmetries, or equations failing the Painleve test. 

Apart from integrability testing, the knowledge of the explicit form of con- 
served densities and higher-order symmetries is useful. For instance, with 
higher-order symmetries of integrable systems, one can build new completely 
integrable systems, or discover connections between integrable equations and 
their group theoretic origin. 

Explicit forms of conserved densities are useful in the numerical solution of 
PDEs or DDEs. In solving DDEs, which may arise from integrable discretiza- 
tions of PDEs, one should check that conserved quantities indeed remain con- 
stant. In particular, the conservation of a positive definite quadratic quantity 
may prevent nonlinear instabilities in the numerical scheme. 

Our integrability package InvariantsSymmetries.m [21] in Mathematica au- 
tomates the tedious computation of closed-form expressions for conserved 
densities and higher-order symmetries for both PDEs and DDEs. Applied 
to systems with parameters, the package determines the conditions on these 
parameters so that a sequence of conserved densities or symmetries exists. 
The software can thus be used to test the integrability of classes of equations 
that model various wave phenomena. Our examples include a vector modified 
KdV equation, the extended Lotka-Volterra and relativitic Toda lattices, and 
the Heisenberg spin model. 

The conserved densities and symmetries presented in this paper were obtained 
with InvariantsSymmetries.m. 
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2 Symbolic Program for the Painleve Test 



2. 1 Purpose 



We focus on PDEs. As originally formulated by Ablowitz et al. [2,3], the 
Painleve conjecture asserts that all similarity reductions of a completely in- 
tegrable PDE should be of Painleve-type; i.e. its solutions should have no 
movable singularities other than "poles" in the complex plane. 

A later version of the Painleve test due to Weiss et al. [52] allows testing of the 
PDE directly, without recourse to the reduction(s) to an ODE. A PDE is said 
to have the Painleve property [1] if its solutions in the complex plane are single- 
valued in the neighborhood of all its movable singularities. In other words, the 
equation must have a solution without any branching around the singular 
points whose positions depend on the initial conditions. For ODEs, it suffices 
to show that the general solution has no worse singularities than movable 
poles, or that no branching occurs around movable essential singularities. 

A three step-algorithm, known as the Painleve test, allows one to verify 
whether or not a given nonlinear system of ODEs or PDEs with (real) polyno- 
mial terms fulfills the necessary conditions for having the Painleve property. 
Such equations are prime candidates for being completely integrable. 

There is a vast amount of literature about the test and its applications to spe- 
cific ODEs and PDEs. Several well-documented surveys [5,9,15,31,32,34,38,40] 
and books [8,10,48] discuss subtleties and pathological cases of the test that 
are far beyond the scope of this article. Other survey papers [15,16,39] deal 
with the many interesting by-products of the Painleve test. For example, they 
show how a truncated Laurent series expansion of the type introduced below, 
allows one to construct Lax pairs, Backhand and Darboux transformations, 
and closed-form particular solutions of PDEs. 



2.2 Algorithm for a Single Equation 



We briefly outline the three steps of the Painleve test for a single PDE, 

F(x,t,u(x,t)) = 0, (1) 
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in two independent variables x and t. Our software can handle the four inde- 
pendent variables (x,y,z,t). Throughout the paper we will use the notations 



du 

Ut "777; U nx 



d 2 u 



(2) 



dx n ' 



Utx — 



dtdx : 



etc. 



In the approach proposed by Weiss, the solution u(x, t), expressed as a Laurent 
series 



should be single-valued in the neighborhood of a non-characteristic, movable 
singular manifold g(x,t), which can be viewed as the surface of the movable 
poles in the complex plane. In (3), u (x,t) 7^ 0, a is a negative integer, and 
Uk(x,t) are analytic functions in a neighborhood of g(x,t). 

Note that for ODEs the singular manifold is g(x,t) = x — xq, where xo is the 
initial value for x. For PDEs, if u(x, t) has simple zeros and g x (x,t) ^ 0, one 
may apply the implicit function theorem near the singularity manifold and 
set g(x, t) = x — h(t), for an arbitrary function h(t) [33,44]. This so-called the 
Kruskal simplification, considerably reduces the length of the calculations. 

The Painleve test proceeds in three steps: 

Step 1: Determine the dominant behavior 

Determine the negative integer a and Uq from the leading order "ansatz" . This 
is done by balancing the minimal power terms after substitution of u oc Uog a 
into the given PDE. There may be several branches for u , and for each the 
next two steps must be performed. 

Step 2: Determine the resonances 

For a selected a and u , calculate the non- negative integer powers r, called 
the resonances, at which arbitrary functions u r enter the expansion. This is 
done by requiring that u r is arbitrary after substitution of u oc uog a + u r g a+r 
into the equation, only retaining its most singular terms. The coefficient u r 
will be arbitrary if its coefficient equals zero. The integer roots of the resulting 
polynomial must be computed. The number of roots, including r = — 1, should 
match the order of the given PDE. The root r = — 1 corresponds to the 
arbitrariness of the manifold g(x,t), 

Step 3: Verify the correct number of free coefficients 



oo 



u(x, t) = g a (x, t) M x , t) g k (x, t), 



(3) 



k=0 
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Verify that the correct number of arbitrary functions u r indeed exists by sub- 
stituting the truncated expansion 



into the PDE, where rmax is the largest resonance. At non-resonance levels, 
determine all u^. At resonance levels, u r should be arbitrary, and since we are 
dealing with a nonlinear equation, a compatibility condition must be verified. 

An equation for which the above steps can be carried out consistently and 
unambiguously, is said to have the Painleve property and is conjectured to be 
completely integrable. This entails that the solution has the necessary number 
of free coefficients u r , and that the compatibility condition at each of these 
resonances is unconditionally satisfied. 

The reader should be warned that the above algorithm does not detect es- 
sential singularities and therefore cannot determine whether or not branching 
occurs about these. So, for an equation to be integrable it is not sufficient that 
it passes the Painleve test. Neither it is necessary. Indeed, there are integrable 
equations, such as the Dym-Kruskal equation, u t = u 3 U3 X , that do not pass the 
Painleve test, yet, by a complicated change of variables can be transformed 
into an integrable equation. 

2.3 Algorithm for Systems 

The generalization of the algorithm to systems of ODEs and PDEs is obvi- 
ous. Yet, it is non-trivial to implement. One of the reasons is that the major 
symbolic packages do not handle inequalities well. 

With respect to systems, our code is based on the above three step-algorithm 
but generalized to systems, as it can be found in [33,44,48]. In these papers 
there is an abundance of worked examples that served as test cases. 

For example, given a system of first-order ODEs, 



rmax 




(4) 




(5) 




one introduces a Laurent series for every dependent variable Ui(x) : 



OG 




(6) 
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The computer program must carefully determine all branches of dominant 
behavior corresponding to various choices of ctj and/or Uq\ For each branch, 
the single-valuedness of the corresponding Laurent expansion must be tested 
(i.e. the resonances must be computed and the compatibility conditions must 
be verified). All the details can be found in [15,33,44]. 

Singularity analysis for PDEs is nontrivial [32] and the Painleve test should 
be applied with extreme care. Notwithstanding, our software automatically 
performs the formal steps of the Painleve test for systems of ODEs and PDEs. 
The examples in section 3.2 illustrate how the code works. Careful analysis 
of the output and drawing conclusions about integrability should be done by 
humans. Some subtleties of the mathematics of the Painleve test of systems 
of PDEs were also dealt with in [6,13,28,49]. 



3 Examples of the Painleve Test 

3.1 Single Equations 

Numerous examples of the Painleve test for ODEs can be found in the re- 
view papers. We turn our attention to PDEs. Using our software package 
painsing. max or painsing.mone can determine the conditions under which the 
equation 

u tx + a(t)u x + 6uu 2x + 6u 2 x + u 4x = 0, (7) 
passes the Painleve test. 

For (7), a = — 2 and u = —2g%. Apart from r = — 1, the roots are r = 4,5, 
and 6. The latter three are resonances. Furthermore, 

Mi = 2g 2x , u 2 = --^(4:g x g 3x - 3gj x + g t g x ), (8) 
vy x 

u 3 = J^(a(t)gl + g 2 x g± x - 4g x g 2x g 3x + 3g 2x - g t g x gi x + gtx&), (9) 

and w 4 and u 5 are indeed arbitrary since the compatibility conditions at res- 
onances r = 4 and r = 5 are satisfied identically. 

The compatibility condition at resonance r = 6 is at + 2a 2 = 0. Ignoring the 
trivial solution, we get a = 2 (t-t ) ' Without loss of generality, we set t = and 
equation (7) becomes the cylindrical KdV equation which is indeed completely 
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integrable [1,4]. Painleve based investigations for integrable PDEs with space 
and time dependent coefficients are given in [1,7,24,27]. 



Our Painleve programs cannot automatically test a class of equations such as 



u t + au 2 u x + bu x u 2x + cuu 3x + u 5x = 0, (10) 



with arbitrary (non-zero and real) parameters a, b and c. The parameters affect 
the lowest coefficient in the Laurent expansion in such as way that the roots 
(r) cannot be computed, and the integrability conditions can no longer be 
tested. 



In (10) there are four cases that are of particular interest: 



(i) a = yqC 2 and b = 2c (Lax equation), 

(ii) a = |c 2 and b = c (Sawada-Kotera equation), 

(iii) a = |c 2 and b — |c (Kaup-Kupershmidt equation), and 

(iv) a = |c 2 and b = 2c (Ito equation). 



In Table 1 we list the results of the Painleve test applied to these cases. For 
the first three equations the compatibility conditions are satisfied at all the 
resonances. These equations pass the test. For the Ito equation the compati- 
bility conditions are only satisfied at some of the resonances. The Ito equation 
fails the test. The first three equations are known to be completely integrable. 
Ito's equation is not completely integrable. 



The two other algorithms presented in this paper can determine the conditions 
(i), (ii) and (iii) that assure the complete integrability of (10). The conserved 
densities and higher-order symmetries of (10) can be found in [18] and [20]. 
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Lax 


Sawada-Kotera 


Kaup— Kupershmidt 


Ito 


(a, 6, c)=(30, 20, 10) 


(a, 6, c) = (5,5,5) 


(a, b, c) = (20,25,10) 


(a,M = (2,6,3) 


a = -2 


a = -2 


a = -2 


a = -2 


Branch 1 


Branch 1 


Branch 1 


Branch 1 


u = -2g 2 x 


mo = -Qg 2 x 


M o = -\gl 


u = -Qg 2 x 


r=-l,2,5,6,8 


r=-l,2,3,6,10 


r=— 1, 3, 5, 6, 7 


r=-l,3,4,6,8 


OK for r > 2 


OK for r > 2 


OK for r > 3 


OK at r = 3 








Not atr = 4,6,8 


Branch 2 


Branch 2 


Branch 2 


Branch 2 


u o = -6^ 


u = -Ylgl 


Mo = -12^ 


«o = -30^ 


r=-3 -1,6,8,10 


r=-2 -1,5,6,12 


r=-7,-l,6,10, 12 


r=-5,-l,6,8,12 


OK for r > 6 


OK for r > 5 


OK for r > 6 


OK at r = 6,8 








Not at r = 12 


Passes Test 


Passes Test 


Passes Test 


Fails Test 



Table 1: Painleve analysis of fifth-order KdV equation 

u t + au 2 u x + bu x u 2x + cuu 3x + u 5x = 



5.^ Simple Systems 

We start with a famous system of ODEs, 

u[ = a(u2 — Ml), «2 = — U1U3 + 6« 1 — U2, u' 3 = U\U<i — CM3, (11) 

where a, 6, and c are positive constants. System (11) was proposed by the 
meteorologist E. N. Lorenz as a simplified model for atmospheric turbulence 
in a vertical air cell beneath a thunderhead [47]. 

Using our code painsys. max or painsys.m, with the series (6) for % — 1, 2, 3, we 
determine the leading orders ol\ — — 1, cx.2 — Qtz — —2. The first coefficients in 
the series are 

«?> = ±2i, 4 2) = T-, 4 3) = --- (12) 
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The roots are r = —1,2,4. Furthermore, 

1 2 
= T -i(2c-3a+ 1), uf ] = ±2i, uf ] = — (3a - c + 1). (13) 
3 3a 

The compatibility conditions at resonances r = 2 and r = 4 are not satisfied. 
At resonance r = 2 we encounter the condition a(c — 2a)(c+ 3a — 1) = 0. So, 
we consider two cases: 

• For c = 2a, the compatibility condition at r = 4 is not satisfied. 

• For c = 1 — 3a, the compatibility condition at r = 4 is satisfied 
provided that a — |. 

We conclude that the Lorenz system (11) passes the Painleve test when a = | 
and c = 0. This special case may not be relevant in the context of meteorology 
since in (11) the parameters a, 6, c in (11) are supposed to be positive. 

To illustrate the generalization of the algorithm to systems of PDEs, consider 
the system 

uu xt - u x u t + v x v t = 0, uv xt - u t v x - u x v t = 0, (14) 

which plays a role in elementary particle physics [49]. For any of the dependent 
variables we introduce a Laurent expansion 

oo oo 

u = g a J2u k g\ v = g ej2v k g k . (15) 

k=0 k=0 

Analysis of the dominant behavior leads to a — f3 — — 1 and Uq + v q = 0. 
Let Vq be arbitrary, then Uq = ±Wq. Next, we look for powers of g at which 
arbitrary functions u r ,v r can enter. Equating the coefficients of g r ~ 4 ,r > 0, 
we obtain the system 

2r ir(r — 1) 

i(r 2 -r + 2) -2(r - 1) 

The determinant of the coefficient matrix vanishes provided r = —1,0,1,2. 
Note that r = confirms that v can be chosen freely. Upon substitution of 
the Laurent series (truncated at level k — 2) and setting the terms in g~ 2 
and g~ 3 equal to zero, one finds that u\ is completely determined, whereas v\ 
and U2 (or i^) are arbitrary. The compatibility conditions are satisfied. The 
Laurent expansions have the required number of arbitrary functions u k - Thus, 
the system (14) passes the Painleve test. 



u r 



(16) 
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Hirota and Satsuma [1] proposed a coupled system of KdV equations, 

u t - 6auu x + 6vv x - au 3x = 0, v t + 3uv x + v 3x = 0, (17) 

where a is a nonzero parameter. System (17) describes interactions of two 
long waves with different dispersion relations. It is known to be completely 
integrable for a — \. This is confirmed by the Painleve test. Indeed, with (15) 
we obtain a — f3 — —2 and r = —2, —1, 3, 4, 6 and 8. Furthermore, uq = —4 
and vo = ±2y / 2a, determine the coefficients u\, i>i, u 2 , v 2 unambiguously. At 
resonances 3 and 4 there is one free function and no condition for a. The coef- 
ficients W5 and t> 5 are unique determined, but at resonance 6, the compatibility 
condition is only satisfied if a = \. For this value, the compatibility condition 
at resonance 8 is also satisfied. 

For the integrable version of the Boussinesq system [45] 

u t + v x + uu x = 0, v t + u 3x + (uv) x = 0, (18) 

with the Laurent expansions in (15), the leading order is a = —1,(3 = —2. 
Careful investigation of the recursion relations linking the Uk and Vk allows one 
to conclude that there are resonances at levels 2, 3 and 4. Finally, the three 
compatibility conditions are seen to hold after lengthy computations. System 
(18) thus passes the test. 

Our codes painsys.max are also applicable to complex equations such as the 
nonlinear Schrodinger (NLS) equation [1], 

iqt - q2x + 2|g| 2 g = 0, (19) 

where q(x,t) is a complex function. First, (19) must be rewritten as a system: 
ut - v 2x + 2v(u 2 + v 2 ) = 0, v t + u 2x - 2u(u 2 + v 2 ) = 0, (20) 



where q = u + iv. For simplicity, we use the Kruskal simplification, setting 
g(x,t) = x — h(t) in (15). Then, the leading order is a = (3 = —1, and 

r = —1, 0, 3 and 4. Furthermore, v = ±v/l — Uq, with uq arbitrary, and 



11 1 

ui = ^v h t , vi = --uoht, u 2 = -j^-u (v h 2 + 2u ,t), 

V2 = -j^(v h 2 + 2uo, t ), u 3 = --^( h tt + 8v v 3 ), 

M 4 = -7^( 2v o h thtt ~ vlu ,tt - u u 2 t - 24u VqV A ), (21) 
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where 1*3 and V4 are arbitrary. Note that at every resonance there is one free 
function in the Laurent expansion. The NLS equation passes the test. 

In searching the literature we found numerous examples of systems for which 
the Painleve test was carried out by hand. Some of the most elaborate exam- 
ples of Painleve analysis involve three- wave interactions [37], the mixmaster 
universe model [12], and chaotic star pulsations [51]. 



4 Symbolic Programs for Conserved Densities and Symmetries 

4-1 The Key Concepts 

The key observation behind our algorithms is that conserved densities and 
higher-order symmetries of a PDE or DDE system abide by the dilation sym- 
metry of that system. We illustrate this for prototypical examples of single 
PDEs and DDEs. 

Dilation Invariance of PDEs. The ubiquitous Korteweg-de Vries (KdV) 
equation from soliton theory [1], 

u t = Quu x + u 3x , (22) 

is invariant under the dilation (scaling) symmetry 

(f , x, u) -> (X'H, X^x, X 2 u), (23) 

where A is an arbitrary parameter. Obviously, u corresponds to two derivatives 
in x, i.e. u ~ d 2 /dx 2 . Similarly, d/dt ~ d 3 /dx 3 . Introducing weights, denoted 
by w, we have w(u) = 2 and w{d/dt) = 3, provided we set w{d/dx) = 1. The 
rank R of a monomial equals the sum of all of its weights. Observe that (22) 
is uniform in rank since all the terms have rank R = 5. 

Conserved Densities of PDEs. For PDEs like (22), the conservation law 

D t p + D,J = (24) 

connects the conserved density p and the associated flux J. As usual, D t and 
D x are total derivatives. With few exceptions, polynomial density-flux pairs 
only depend on u, u^,, etc., and not explicitly on t and x. 

The first three (of infinitely many) independent conservation laws for (22) are 
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D t (u)-D x {3u 2 + u 2x ) = 0, D t (u 2 ) - B x (Au 3 - u\ + 2uu 2x ) = 0, (25) 

19 1 

D t (u 3 - -u 2 x ) - ^x{-u A - Quu 2 x + 3u 2 u 2x + -u 2 2x - u x u 3x ) = 0. (26) 

The densities p = u,u 2 ,u 3 — \u 2 £ have ranks 2,4 and 6, respectively. The 
associated fluxes have ranks 4, 6 and 8. The terms in the conservation laws 
have ranks 5, 7, and 9. 

Integration of both terms in the conservation law with respect to x yields 

+oo 

P = J p dx = constant in time, (27) 

— oo 

provided J vanishes at infinity. P is the true conserved quantity. The first two 
conservation laws correspond to conservation of momentum and energy. For 
ODEs, the quantities P are called constants of motion. 

Symmetries of PDEs. As summarized in Table 2, G(x, t, u, u x , u 2a; , ...) is 
a symmetry of a PDE system if and only if it leaves it invariant for the change 
u — > u + eG within order e. Hence, D t (u + eG) = F(u + eG) must hold up to 
order e. Thus, G must satisfy the linearized equation D t G = F'(u)[G], where 
F' is the Frechet derivative: F'(u)[G] = £ F(u + eG)| e=0 . 





Continuous Case (PDEs) 


Semi-discrete Case (DDEs) 


System 


u t = F(u, u x ,u 2x ,...) 


u n = F(..., U n _i, U n , U n+ i, ...) 


Conservation Law 


D tP +D x J = 


P n + Jn+1 ~ Jn = 


Symmetry 


D t G = F'(u)[G] 

= |F(u + eG)| e=0 


D t G = F'(u n )[G] 

= |F(u n + e G)| e=0 



Table 2: Conservation Laws and Symmetries 
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Two nontrivial higher-order symmetries [42] of (22) are 



G {1) =30u 2 u x + 20u x u 2x + 10uu 3x + u 5x , (28) 
G ,(2) =140u 3 u a; + 70u 3 x + 280uu x u 2x + 70u 2 u 3x + 70u 2x u 3x + A2u x u Ax 

+lAuu 5x + u 7x . (29) 

The recursion operator [42, p. 312], 

$ = D 2 + Au + 2u x D-\ (30) 

connects the above symmetries, $ = , and also the lower order sym- 
metries as shown in [20]. 

Note that the recursion operator (30) is also uniform in rank with R = 2 since 
w{D~ 1 ) — — 1. Currently, we are working on an algorithm and symbolic code 
to compute recursion operators based on the knowledge of a few higher-order 
symmetries and conserved densities. 

Higher-order symmetries lead to new integrable PDEs. Indeed, the evolution 
equations u t = G^ and u t = G^ are the well-known fifth and seventh-order 
equations in the completely integrable KdV hierarchy [1]. 

Dilation Invariance of DDEs. We now turn to the semi-discrete case. As 
prototype, consider the Volterra lattice [1], 

u n = u n (u n+ i - u n -i). (31) 

which is one of the discretizations of (22). Note that (31) is invariant under 

{t,u n )^(\-H,\u n ). (32) 

So, u n ~ d/dt, or w(u n ) = 1 if we set u>(d/dt) = 1. Every term in (31) has 
rank R = 2, thus (31) is uniform in rank. 

Conserved Densities of DDEs. For DDEs like (31), the conservation law 

Pn + Jn+1 ~ J n = (33) 



connects the conserved density p n and the associated flux J n . For (31) the 
first two conservation laws (of ranks 2 and 3) are 



d 

— {u n ) + u n u n _ x - u n+1 u n = 0, 



(34) 
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^2 Un + UnUn+1 ^ + Un ~ lU n + «n-lMTi« n +l - U n U n+1 - U n U n+1 U n+2 = 0. (35) 

The densities p n = u n and p n = u n {\u n + u n+1 ) have ranks 1 and 2. Their 
fluxes J n = —u n u n ^i and J n = —u n -\U 2 n — u n -iU n u n+ \ have ranks 2 and 3. 

Symmetries of DDEs. For DDEs of type (31), G(..., u„_i, u„, u n+1 , ...) is 
a symmetry if and only if the infinitesimal transformation u„ — > u„ + eG 
leaves the DDE invariant within order e. Consequently, G must satisfy ^ = 
F'(u„)[G], where F' is the Frechet derivative, F'(u„)[G] = £F(u n + eG)| £=0 . 
See Table 2. 

The first nontrivial higher-order symmetry of (31) is 

G = U n U n+ i (u n + U n+ i + U n+2 ) - U n -iU n (w n _ 2 + Un-l + Un), (36) 

and, similar to the continuous case, an integrable lattice u n = G follows. 

Both (22) and (31) have infinitely many polynomial conserved densities [18,22] 
and symmetries [20]. 

Remarks: 

(i) For scaling invariant systems like (22) and (31), it suffices to consider the 
dilation symmetry on the space of independent and dependent variables. 

(ii) For systems that are inhomogeneous under a suitable scaling symmetry 
we use the following trick: We introduce one (or more) auxiliary parameter(s) 
with appropriate scaling. In other words, we extend the action of the dilation 
symmetry to the space of independent and dependent variables, including the 
parameters. These extra parameters should be viewed as additional dependent 
variables, with the caveat that their derivatives are zero. With this trick we can 
apply our algorithms to a larger class of polynomial PDE and DDE systems. 
Examples can be found in [18-20,22]. 

Scaling (dilation) invariance, which is a special Lie-point symmetry, is common 
to many integrable nonlinear PDEs and DDEs. In the next two sections we 
show how the scaling invariance can be explicitly used to compute polynomial 
conserved densities and higher-order symmetries of PDEs and DDEs. 

4-2 Algorithms for the Computation of Conserved Densities 

PDE Case. Conserved densities of PDEs can be computed as follows: 
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• Require that each equation in the system of PDEs is uniform in rank. Solve 
the resulting linear system to determine the weights of the dependent variables. 
For instance for (22), solve w(u) + w(d/dt) = 2w(u) + 1 = w(u) + 3, to get 
w(u) = 2 and w(d/dt) = 3. 

• Select the rank R of p, say, R = 6. Make a linear combination of all the 
monomials in the components of u and their x-derivatives that have rank R. 
Remove all monomials that are total ^-derivatives (like u^). Remove equiv- 
alent monomials, that is, those that only differ by a total x-derivative. For 
example, uu<i x and u x are equivalent since UU2 X = \{u 2 )2 X — u x - For (22), one 
gets p = c\u 3 + C2U 2 X of rank R = 6. 

• Substitute p into the conservation law (24). Use the PDE system to eliminate 
all t-derivatives, and require that the resulting expression E is a total x- 
derivative. Apply the Euler operator (see [18] for the general form), 

£ « - 1 ~ D '^J + D -^J + " ' + (37) 

to E to avoid integration by parts. If any terms remain, they must vanish 
identically. This yields a linear system for the constants q. Solve the system. 
For (22), one gets c l = l,c 2 = -1/2. 

See [18] for the complete algorithm and its implementation. 

DDE Case. The computation of conserved densities proceeds as follows: 

• Compute the weights in the same way as for PDEs. For (31), one gets 
w(u n ) = 1 by solving w(u n ) + 1 = 2w(u n ). Note that u>(d/dt) = 1 and weights 
are independent of n. 

• Determine all monomials of rank R in the components of u„ and their 
t-derivatives. Use the DDE to replace all the t-derivatives. Monomials are 
equivalent if they belong to the same equivalence class of shifted monomials. 
For example, u n -iV n+ i,u n+ 2V n+ 4 and -u n _3f„_i are equivalent. Keep only the 
main representatives (centered at n) of the various classes. 

• Combine these representatives linearly with coefficients q, and substitute 
the form of p n into the conservation law p n — J n — J n+1 . 

• Remove all t-derivatives and pattern-match the resulting expression with 
Jn ~ Jn+i- To do so use the following equivalence criterion: if two monomials 
mi and vri2 are equivalent, m\ = m.2, then mi = m 2 + [M n — M n+ i] for some 
polynomial M n that depends on u n and its shifts. For example, u n -2U n = 

Un-lUn+l Since U n _ 2 U n = W„-lW n+ i + [« ri -2«n-«n-lMn+l] = Mn-l«n+l + [M n - 
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M n+ i], with M n = u n -2U n . Set the non- matching part equal to zero, and solve 
the linear system for the q. Determine J n from the pattern J n — J n+ ±. For 
(31), the first three (of infinitely many) densities p n are listed in Table 3. 





KdV Equation 


Volterra Lattice 


Equation 


u t = 6uu x + u 3x 


u n = u n (u n+1 - U n -i) 


Densities 


p = u, p = u 2 
p = u 3 - \u\ 


Pn — U n , p n = U n (^U n + U n+ i) 

Pn = \ul+ U n U n+1 (u n + U n+1 +U n+2 ) 


Symmetries 


G = u x , G = 6uu x + u 3x 
G = 30u 2 u x + 20u x U2 X 
+10uu 3x + u 5x 


G = u n u n+1 (u n + u n+ i + u n+2 ) 

-U n -iU n {u n - 2 + U n -i + U n ) 



Table 3: Prototypical Examples 



Details about this algorithm and its implementation are in [19,22]. See [21] for 
an integrated Mathematica Package that computes conserved densities (and 
also symmetries) of PDEs and DDEs. 

4-3 Algorithm for the Computation of Symmetries 

PDE Case. Higher-order (or generalized symmetries) of PDEs can be com- 
puted as follows: 

• Determine the weights of the dependent variables in the system. 

• Select the rank R of the symmetry. Make a linear combination of all the 
monomials involving u and its ^-derivatives of rank R. For example, for (22), 
G — Ci u 2 u x +C2 u x u 2x +c 3 uusx+Ci Ur JX is the form of the generalized symmetry 
of rank R = 7. In contrast to the computation of conserved densities, no terms 
are removed here. 

• Compute D(G. Use the PDE system to remove all t-derivatives. Equate the 
result to the Frechet derivative F'(u)[G]. Treat the different monomial terms 
in u and its ^-derivatives as independent, to get the linear system for q. Solve 
that system. For (22), one obtains 

G = 30u 2 u x + 20u x u 2x + I0uu 3x + u 5x . (38) 

The symmetries of the Lax family of rank 3, 5, and 7 are listed in Table 3. 
They are the first three of infinitely many. See [20] for the details about the 
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algorithm. 
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DDE Case. Higher-order symmetries of DDEs can be computed as follows: 

• First determine the weights of the variables in the DDE the same way as 
for conserved densities. 

• Determine all monomials of rank R in the components of u„ and their 
t-derivatives. Use the DDE to replace all the t-derivatives. Make a linear com- 
bination of the resulting monomials with coefficients q. 

• Compute D t G and remove all u n -i, u n , Un+i, etc Equate the resulting ex- 
pression to the Frechet derivative F'(u n )[G] and solve the system for the q, 
treating the monomials in u n and its shifts as independent. 

For (31), the symmetry G of rank R = 3 is listed in Table 3. There are 
infinitely many symmetries, all with different ranks. 

See [20] for the complete algorithm and its implementation in Mathematica, 
and [21] for an integrated Mathematica Package that computes symmetries of 
PDEs and DDEs. 

Notes: 

(i) A slight modification of these methods allows one to find conserved densities 
and symmetries of PDEs that explicitly depend on t and x. See the first 
example in the next section. 

(ii) Applied to systems with free parameters, the linear system for the q 
will depend on these parameters. A careful analysis of the eliminant leads 
to conditions on these parameters so that a sequence of conserved densities 
or symmetries exists. Details about this type of analysis and its computer 
implementation can be found in [18]. 



5 Examples of Densities and Symmetries 

5.1 Vector Modified KdV Equation 

In [50], Verheest investigated the integrability of a vector form of the modified 
KdV equation (vmKdV), 




(39) 



or component-wise for B = (w, v), 
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u t + 3u 2 u x + v 2 u x + 2uvv x + us x = 0, 

v t + 3v 2 v x + u 2 v x + 2uvu x + v 3x = 0. (40) 
With our software InvariantsSymmetries.m [21] we computed 

pi = u, p 2 = v, p 3 = u 2 + v 2 , p 4 = ^(u 2 + v 2 ) 2 - (v 2 +v 2 x ), (41) 
P5 = \x{u 2 + v 2 ) - l -t{u 2 + v 2 ) 2 + t(u 2 x + v 2 x ). (42) 

Note that the latter density depends explicitly on x and t. Verheest [50] has 
shown that (40) is non-integrable for it lacks a bi-Hamiltonian structure and 
recursion operator. We were unable to find additional polynomial conserved 
densities. Polynomial higher-order symmetries for (40) do not appear to exist. 

5.2 Heisenberg Spin Model 

The continuous Heisenberg spin system [14] or Landau-Lifshitz equation, 

S t = S x AS + S x DS, (43) 

models a continuous anisotropic Heisenberg ferromagnet. It is considered a 
universal integrable system since various known integrable PDEs, such as the 
NLS and sine-Gordon equations, can be derived from it. In (43), S = [u, v, w] T 
with real components, A = V 2 is the Laplacian, D is a diagonal matrix, and 
x is the standard cross product of vectors. 

Split into components, (43) reads 

u t = vw 2x - wv 2x + {(3 - a)vw, 
v t = wu 2x - uw 2x + (1 - f3)uw, 

w t = uv 2x - vu 2x + (a - l)uv . (44) 

In Table 4 we list the conserved densities for three typical cases; other cases 
are similar. 

Note that for all the cases we considered 

p = M 2 +u 2 +u; 2 = ||g||2 (45) 

is constant in time (since J = 0). Hence, all even powers of ||S|| are also 
conserved densities, but they are dependent of (45). 
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D = diag(l, a, j3) 


D = diag(l, a, 0) 
a ^ 


D = diag(0, 0,0) 


p — u if a — (5 

p = vifp = l 

p = w if a = 1 

p = u 2 + v 2 + w 2 

p= (1 - a )v 2 + (1 - (5)w 2 

+ u l + v l + w l 


p — w if a — 1 

p = M 2 + f 2 + w 2 
p = (1 — a)-u 2 + w 2 


p — u 

p = V 

p = w 

p = u 2 + v 2 + w 2 
P = u 2 x + v 2 x + w 2 x 



Table 4: Conserved Densities for the Heisenberg Spin Model 

Furthermore, the sum of two conserved densities is a conserved density. Hence, 
after adding (45), 



p = (a- l)v 2 + (P- l)w 2 - (u 2 x + v 2 x + w 2 x ) 
can be replaced by 

p — u 2 + oiv 2 4" (3w 2 - (w 2 +v 2 x + w 2 x ). 



(46) 



(47) 



Note that u 2 = D x {uu x ) —uu 2x and recall that densities are equivalent if they 
only differ by a total re— derivative. So, (47) is equivalent with 



p = u 2 + av 2 + (3w 2 + uu 2x + vv 2x + w 2x , 



(48) 



which can be compactly written as p = S-AS+S-ZJS, where D = diag(l, a, (3). 
Consequently, the Hamiltonian of (43) 

H = - l - J S • AS + S • DS dx (49) 
is constant in time. The dot (•) refers to the standard inner product of vectors. 
5.3 Extended Lotka-Volterra and Relativistic Toda Lattices 



Itoh [29] studied this extended version of the Lotka-Volterra equation (31): 



fe-i 



(50) 



r=l 
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For k = 2, (50) is (31), for which three conserved densities and one symmetry 
are listed in Table 3. In [19], we gave two additional densities and in [20] we 
listed two more symmetries. 

For (50), we computed 5 densities and 2 higher-order symmetries for k = 3 
through k = 5. Here is a partial list of the results: 

For k = 3 : 

Pi = u n , P2 = ^u 2 n + u n (u n+1 + u n+2 ), (51) 

1 



P3 = + u 2 n {u n+ i + u n+2 ) + u n {u n+ i + u n+2 ) 2 

+u n (u n+1 u n+3 + u n+2 u n+3 + u n+2 u n+4 ), (52) 

G = U 2 n {u n+ i + U n+2 - M n _ 2 - ttn-l) + U n [(u n+1 + U n+2 ) 2 

-(u n -2 + Mn-l) 2 ] + U n [u n+1 U n+3 + U n+2 U n+3 + U n+2 U n+4 
-{Un-4U n -2 + U n _ 3 U n _ 2 + U n _ 3 U n _i)}. (53) 

For k = 4 : 

Pi = u n , p 2 = \i 2 n + u n (u n+1 + u n+2 + u n+3 ), (54) 

X 

p 3 = -ul + u 2 n (u n+ i + u n+2 + u n+3 ) + u n (u n+ i + u n+2 + u n+3 ) 2 

+U n {u n+ iU n+ A + U n+2 U n+ 4 + U n+3 U n+ 4 + U n+2 U n+ 5 

+u n+3 u n+5 + u n+3 u n+6 ), (55) 

G = W n [ttri+l'Un+4 + u n+2U n +4: + tt n+ 3U n+ 4 + u n+2 u n+5 

+U n+3 U n+5 + U n+3 U n+6 - (-U n _ 6 -U„_ 3 + U n sU n - 3 + It n _4'U n _3 

+U n _5U n _2 - -U„-4Wn-2 + W n -4Wn-l)] + M„[(m„+1 + U n+2 + M n+3 ) 2 

-M„(M n _3 + M„_ 2 + Mn-l) 2 ] + U 2 n [u n+ i + W„ +2 + W«+3 

-(« n _ 3 + 14 n _2 + «„-i)]. (56) 

Our last example involves the integrable relativistic Toda lattice [41]: 

tt n = U n (V n +1 - V n + U n+1 - U n -i), V n = V n (u n - U„_i) . (57) 



We computed the densities of rank 1 through 5. The first three are 



Pi = u n + v n , p 2 = \(u 2 n + w 2 ) + u n (u n+1 + v n + v n+1 ), (58) 

P3 = |(«n + U n) + M n( M n+l + + «n+l) + «n[(«n+l + *Vfl) 2 
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+U n+1 U n+2 + Un+lVn + U n+ iV n+2 + V n V n+1 + V^}. (59) 

We computed the symmetries for ranks (2, 2) through (4, 4). The first two are: 




U n (u n+1 - M„_i + V n+1 - V n ), G ( 2 ] =V n (u n - Un-i), 

U 2 n (u n+1 - U n _ 1 + V n+l - V n )+ U n [(u n+l + V n+ if - (u n _! - 
+U n+1 (u n+2 + V n+2 ) - U n -i(u n - 2 + V n -i)], 

V 2 n {u n - U n -i) + V n (u 2 n - U n -iU n -2 ~ «n-l + 
-Un^Vn-! + U n V n+1 ). 




(60) 



(61) 




(62) 



Conserved densities and symmetries of other relativistic lattices are in [20,22]. 



6 Conclusions 



We presented three methods to test the integrability of differential equations 
and difference-differential equations. One of these methods is the Painleve 
test, which is applicable to polynomial systems of ODEs and PDEs. 

The two other methods are based on the principle of dilation invariance. Thus 
far, they can only be applied to polynomial systems of evolution equations. 
As shown, it is easy to adapt these methods to the DDE case. 

Although restricted to polynomial equations, the techniques presented in this 
paper are algorithmic and have the advantage that they are fairly easy to 
implement in symbolic code. 

Applied to systems with parameters, the codes allow one to determine the 
conditions on the parameters so that the systems pass the Painleve test, or 
have a sequence of conserved densities or higher-order symmetries. Given a 
class of equations, the software can thus be used to pick out the candidates 
for complete integrability. 

Currently, we are extending our algorithms to the symbolic computation of 
recursion operators of evolution equations. In the future we will investigate 
generalizations of our methods to PDEs and DDEs in multiple space dimen- 
sions. The potential use of Lie-point symmetries other than dilation (scaling) 
symmetries will also be studied. 
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